!-------------------------------------------------------------------------------------
subroutine io_read_bgk(fbgk)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only  :  runname, Vmax, pi

  implicit none

  real, intent(inout)     :: fbgk(jv, jz)
  real                    :: dfbgk(jv, jz)

  real                    :: f0
  integer                 :: iv

  open(90, file= trim(runname)//'.bgk', &
       form='unformatted', access='stream')
  read(90) dfbgk
  close(90)

  fbgk = fbgk + dfbgk
  
end subroutine io_read_bgk


!-------------------------------------------------------------------------------------

subroutine set_ic( dfu, uxMat, Fbgk, it)    ! outputs

  use dimensions
  use mpi_wrap
  use fft_3d
 
  use parameters, only : ux0, pi, restart, runname, kz, ModAmp

  implicit none


  real,    intent(out) :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(out) :: uxMat(jz, jy, jx, jDir)
  real,    intent(in)  :: Fbgk(jv,jz)
  integer, intent(out) :: it

  integer              :: ix, iz, iv, iDir, filenum
  character*7 :: file_ext

  real        :: epsilon
  complex     :: noise(jz,jy,jx)

!-- read restart or create new --

  if (restart) then

       write(msg,*) char(10), ' Reading restart', char(10)
       call mpi_wrap_msg

       open(93, file= trim(runname)//'.tag')
       read(93, *) it, filenum 
       close(93)


       file_ext = '0000000'
       write(file_ext,"(i4.4'.'i2.2)") filenum, myid

       open(91, file= trim(runname)//'.dfu.'//file_ext, &
            form='unformatted', access='stream')
       read(91) dfu
       close(91)

       if (jx .gt. 1) then
          open(92, file= trim(runname)//'.ux.'//file_ext, &
               form='unformatted', access='stream')
          read(92) uxMat
          close(92)
       end if

  else
 
       write(msg,*) char(10), ' Creating initial data', char(10)
       call mpi_wrap_msg

       it = 0
       dfu=0.0
       uxMat=0.0 

       !call random_number(dfu)
       !dfu(:,:,:,:,:) = (dfu(:,:,:,:,:)*2 - 1)*ModAmp
       !
       !do iDir=1,jDir
       !   do iv=1,jv
       !      noise=dfu(iv,:,:,:,idir);
       !      call fft_forward_z(noise);
       !      noise(1,:,:)=0;
       !      call fft_backward_z(noise);
       !      dfu(iv,:,:,:,idir)=noise;
       !   enddo
       !enddo

       !-- multi-mode
       !
       !call random_number(uxMat)
       !do iDir = 1, jDir
       !   uxMat(:,:,:,iDir) = (uxMat(:,:,:,iDir)*2 - 1)*ModAmp + ux0(iDir)
       !enddo


       !-- single-mode
       
       noise = 0;
       noise(2,1,2) = ModAmp;
       
       call fft_backward_x(noise);
       call fft_backward_y(noise);
       call fft_backward_z(noise);
       
       do iDir = 1, jDir
          uxMat(:,:,:,iDir) = (1 + noise(:,:,:))*ux0(iDir)
       enddo


       !-- single kx,kz mode

       !do iz = 1,jz
       !   do iv = 1,jv
       !      dfu(iv,iz,:,1,:) =   0.01 * cos(kz*zvec(iz)) * Fbgk(iv,iz)
       !      dfu(iv,iz,:,3,:) =  -0.01 * cos(kz*zvec(iz)) * Fbgk(iv,iz)
       !   end do
       !end do

  end if



end subroutine set_ic

!-------------------------------------------------------------------------------------
subroutine write_bin(dfu, uxMat, Ex2D, Ez2D, it, filenum)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only  :  runname
  use mpi_wrap

  implicit none

  real,    intent(in)    :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)    :: uxMat(jz, jy, jx, jDir)
  real,    intent(in)    :: Ex2D(jz,jy,jx)
  real,    intent(in)    :: Ez2D(jz,jy,jx)
  integer, intent(in)    :: it
  integer, intent(in)    :: filenum


  integer     :: flag
  character*7 :: file_ext

  file_ext = '0000000'
  write(file_ext,"(i4.4'.'i2.2)") filenum, myid


  open(91, file= trim(runname)//'.dfu.'//file_ext, &
         form='unformatted', access='stream')
  write(91) dfu
  close(91)

  if (jx .gt. 1) then
     open(92, file= trim(runname)//'.ux.'//file_ext, &
          form='unformatted', access='stream')
     write(92) uxMat
     close(92)
  end if


  if (myid .eq. 0) then

!    open(92, file= trim(runname)//'.Ex.'//file_ext, &
!          form='unformatted', access='stream')
!     write(92) Ex2D
!     close(92)

!     open(92, file= trim(runname)//'.Ez.'//file_ext, &
!          form='unformatted', access='stream')
!     write(92) Ez2D
!     close(92)


     open(93, file= trim(runname)//'.tag')
     write(93,*) it, filenum 
     close(93)

  end if

 end subroutine write_bin


!-------------------------------------------------------------------------------------
subroutine diagnostics( dfu, Ex2D, Ez2D, uxMat, t)            ! inputs
!-------------------------------------------------------------------------------------

  use dimensions
  use mpi_wrap
  use fft_3d

  use parameters, only : runname, restart, ux0

  implicit none

  !-- global variables --

  real, intent(in)       :: dfu(jv, jz, jy, jx, jDir)
  real, intent(in)       :: Ex2D(jz,jy,jx)
  real, intent(in)       :: Ez2D(jz,jy,jx)
  real, intent(in)       :: uxMat(jz, jy, jx, jDir)
  real, intent(in)       :: t


  !-- local variables --

  complex    :: Ek(jz,jy,jx)
  real       :: Eksq(jz,jx)

  real       :: KE, PotEx, potEz
  real       :: funiform(jv), funiform1(jv)
  real       :: uxsq(jDir), uxsq1(jDir)

  integer    :: ix, iy, iz, iv, iDir

  logical, save          :: first_call = .true.
  integer, save          :: filenum = 0

  character*4            :: file_ext


  !-----------------------------------------------------------------------------------
       
       
      do iv=1,jv
         funiform1(iv) = sum( dfu(iv,:,:,:,:) )
      end do

      do iDir=1,jDir
         uxsq1(iDir) = sum ( ( uxMat(:,:,:,iDir) + ux0(iDir) )**2 )
      end do


      call mpi_reduce(funiform1, funiform, jv, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)

      call mpi_reduce(uxsq1, uxsq, jDir, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)

      if (myid .ne. 0) return
     
      funiform = funiform /(jx*jy*jz)
      uxsq     = uxsq /(jx*jy*jz)


      PotEx = 0.5*sum(Ex2D(:,:,:)**2) / (jx*jy*jz)
      PotEz = 0.5*sum(Ez2D(:,:,:)**2) / (jx*jy*jz)

      KE =  0.5 * sum( funiform(:) * vvec(:)**2 ) * dv

      !-- individual modes of Ex --

      Ek = Ex2D

      call fft_forward_z(Ek)
      call fft_forward_x(Ek)
      call fft_forward_y(Ek)

      Eksq = abs(Ek(:,1,:))
      Eksq = Eksq**2


      !------------------------------
       
      open(3, file=trim(runname)//'.dat', position='append')

      if (first_call .and. .not. restart) then
  
         write(3, "('% 1.t  2.KineticEnergy  3.Ex^2energy  4.Ez^2energy  5.Ux^2average')")
         write(3,*)

      end if

      write(3,501) t, KE, potEx, potEz, sum(uxsq)
 
      close(3)


      !-- individual modes

      filenum = filenum+1
 
!     open(5, file= trim(runname)//'.Ek', &
!          form='unformatted', access='direct',recl=jx*jz*8)
!     write(5, rec=filenum ) Eksq
!     close(5)
  

      !------------------------------

      first_call = .false.

    501   format(  6400(e12.6,2x)  )


end subroutine diagnostics


!-------------------------------------------------------------------------------------
